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ABSTRACT 


In dynamic fracture problem role of material inertia becomes significant. The formulation 
and analytical solution of dynamic fracture problem are of very complex nature. 
Experimental and numerical methods like finite element method are used extensively for 
solving dynamic fracture problem. Existing methods for simulating crack phenomenon in 
dynamic fracture problems have certain drawbacks such as , they are either 
computationally very intensive or give oscillations in the solution. 

In this work, a model is proposed for simulating dynamic crack propagation 
phenomenon in finite element method. Gradual propagation of crack from one node to 
the next node is simulated by adding a 1-D elastic element at the crack tip node, whose 
stiffness value is reduced gradually to zero as the crack propagate in the element. Shape 
function for elemental mass matrix is modified so that the mass of the element 
surrounding the crack tip also varies as a function of crack length within the element, 
starting from zero to the value corresponding to regular shape function. 

The finite element model is first validated with hypothetical data for different crack 
speeds. Results give smooth variation of energy release rate curve for various crack 
velocities. This method is also applied to determine the dynamic Gi for the experimental 
data on glass fabric/epoxy composites. Results show that the energy release rate gives 
fairly smooth variation and stable solution 



List of Figures 


2.1 Crack Opening scheme in force release model 16 

2.2 Discontinuously moving singular element 18 

3.1 Crack opening pattern in model 24 

3.2 F Vs a curve for different C 28 

3.3 Representation of an element at the crack tip 30 

3.4 Shape function Nj for a=l 30 

3.5 Variation of shape function Ni Vs crack length(r] = -1) 32 

3.6 Shape function Ni for a=0.5 32 

3.7 Elemental mass matrix for a=0.5 and a=l .0 34 

4.1(a) DCB specimen (Full) 45 

4.1(b) DCB specimen (One half) 46 

4.2 Input force pulse 46 

4.3 Energy release rate variation for quasi-static crack propagation 47 

4.4 Energy release rate variation for crack velocity 1500m/s and q-s C value 48 

4.5 Effect of change in C on effective stiffness value at the crack tip 49 

4.6 Effect of change in C on energy release rate 50 

4.7 Effect of averaging f on energy release rate 51 

4.8 Energy release rate after shape function modification of elemental 52 

mass matrix 

4.9 Energy release rate for crack velocity of 1250m/s 53 

4.10 Effective stiffness at crack tip Vs ‘a’ without shape function modification 54 

4.1 1 Variation of effective mass at crack tip Vs ‘a’ 55 

4.12 Variation of effective stiffness at crack tip Vs ‘a’ 56 

4.13 Energy release rate for different time steps 57 

4.14 Resultant force using Newmark’s method 58 

4.15 Resultant force using Q -method 58 

4.16 Resultant force after smoothening by Bezier curve 59 

4.17 Displacement of the cantilever end Expt 1 60 

iii 



4.18 Displacement of the cantilever end Expt 2 60 

4.19 Displacement of the cantilever end Expt 3 61 

4.20 Energy release rate for Expt 1 (present study) 62 

4.21 Energy release rate for Expt 2 (present study) 63 

4.22 Energy release rate for Expt 3 (present study) 64 

4.23 Energy release rate for Expt 1 by Force release model [27] 65 

4.24 Energy release rate for Expt 2 by Force release model [27] 66 

4.25 Energy release rate for Expt 3 by Force release model [27] 67 



List of Tables 


4.1 Details of specimen 43 

4.2 Initiation and propagation toughness with ‘Force release model’ and ‘present’. 44 
model 


v 



LIST OF SYMBOLS 


a 

a' 

A 

B 

[B] 

c 

C,Ci,C 2 

[D] 

E 


F, 


HB 


G 


J prop 


K n 


k d 

M 


L 

[M] 


Ne 

M 


Non dimensional crack length 

Crack length 

Crack surface area 

Thickness of the specimen 

Elastic strain displacement matrix 

Dilatation wave velocity 

Constants in stiffness release formulation 

Elastic constitutive relation matrix 

Young’s modulus of elasticity 

Holding back force at crack tip 

Global load vector 

Energy release rate 

Initiation toughness 

Propagation toughness 

Spring stiffness at the crack tip 

Static stiffness at the crack tip 
Damping coefficient 
Global stiffness matrix 
Effective global stiffness matrix 
Length of specimen 
Global mass matrix 
Number of elements 
Matrix of interpolation function 


VI 



{ 4 } 

fel 

At 

T 


Tin 


U 


0 


u 

{u} 

V c 

w 

rr ext 

w 

Z,T1 


n 

p 

V 

{&?} 

{ 5s } 

{“■} 


Velocity component 
Acceleration component 
Time step 
Kinetic energy 
Initiation toughness 

Displacement of node when crack tip reaches the next node and element 

opens up completely. 

Strain energy in body 
Displacement vector 
Crack velocity 

Work done due to external forces 

Width of specimen 
Natural coordinates 

Total potential energy 
Density of material 
Poisson ratio 
Virtual displacements 
Virtual strains 
Stress component 



CHAPTER 1 


INTRODUCTION 


1.1 Introduction 

Despite the development in many fields of fracture mechanics and its numerous 
applications, formulation and analytical solution of dynamic problem of this theory 
remained unexplored, until recently, on account of their extremely intricate nature. A 
large number of problems still remain unsolved in the description of this process on 
microscopic and macroscopic levels. 

Time dependent applied load may arise by vibration of structure or imbalance 
in the machine, wave impact, seismic distribution etc. These kinds of problems are 
often solved by neglecting time dependence and are solved under static or quasi-static 
conditions. Quasi-static mechanics of brittle fracture render only first approximation 
to the description of fracture and can simply indicate whether or not catastrophic 
growth of crack sets in. 

Structural members that are subjected to cyclic loading such as vibrations, 
microscopic imperfections in the material give rise to microscopic cracks, which in 
turn grow into crack of significant size over a passage of time. At this stage further 
growth of crack in a stable or unstable manner depends not only on the material but 
also on the nature of loading i.e. static or dynamic. An understanding of the 
mechanics of dynamic fracture is necessary for developing sound design 
methodologies. 
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1.2 Mechanics of Dynamic Fracture 


Subject of dynamic fracture can be broadly described as that of mechanics of solids, 
containing stationary or propagating cracks, wherein effect of material inertia and 
stress wave interaction play significant role. 

Inertia effect can arise, either from rapidly applied loading on the cracked 
body or from rapid crack propagation. In case of rapid loading, influence of load is 
transferred to the crack by means of stress wave through the material. In case of rapid 
crack propagation material particles on the opposite crack faces displace with respect 
to each other. 

As per the guidelines given by Freund [1], in former case, inertia effect will be 
significant if the characteristic time (i.e. maximum load by rate of load increase) is 
larger than time required for stress wave to travel at the characteristic wave speed of 
the material over a representative length of the body (crack length or distance from 
crack edge to the loaded boundary). 

Fundamental framework of dynamic fracture mechanics relies primarily on solution 
for dynamic behaviour of solid containing crack. These solutions are characterised by 

i) Stress wave integration 

ii) Need to take into account kinetic energy in the global energy balance of the 
fracturing body 

iii) Inertia effect of the material 

One logical way of classifying dynamic fracture problem is as follows (Nishioka and 
Atluri[2] ) 

i) Solids containing stationary cracks subjected to dynamic loading 

ii) Solids containing dynamically propagating cracks under quasi-static loading 

iii) Solids containing dynamically propagating cracks under dynamic loading 
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Analytical solutions to the dynamic fracture problem for some cases are 
available. However, these solutions are limited to cases of simple loading and of 
unbounded plane bodies. However, the interaction of stress waves emanating from 
the crack tip or those reflected from the boundaries make closed form solution of 
dynamic problem intractable and it becomes mandatory to analyse such problem 
using computational methods. 

Even laboratory evaluation of dynamic fracture toughness is based on the 
measurement of instantaneous dynamic stress field close to the propagating crack tip. 
Such direct measurements are difficult and require sophisticated instruments and 
accurate measurements. To overcome these difficulties, hybrid experimental and 
numerical methods are often preferred. This is one of the reasons for the advancement 
of the state of science of dynamic fracture mechanics relies heavily on the 
simultaneous advancement in computational methods. 


1.3 Literature Survey 

Though considerable amount of work has been done to study the fracture phenomena 
through numerical methods, only rectangular double cantilever beam (DCB) 
specimens were given special interest. Owen and Shantaram [3] started the case of 
finite element method to study dynamic crack growth. They studied DCB specimen 
and pipeline problem under transient loading. The crack was advanced from one node 
to the next node and when stress at the gauss point nearest to the crack tip exceeded a 
certain value, damping coefficient was made zero at the released node. The crack 
propagation history was simulated for the given loading conditions. 

Nishioka and Atluri [4] investigated the crack propagation and arrest in high 
strength steel DCB specimen using moving singular dynamic finite element 
procedure. An edge crack in a rectangular DCB specimen was propagated by 


3 



procedure. An edge crack in a rectangular DCB specimen was propagated by 
inserting a wedge and the results were compared with experimental data. In another 
work, Nishioka and Atluri [5] presented the results of generation and prediction 
studies of dynamic crack propagation in plane stress and plane strain cases. The 
studies were conducted by using FEM, taking into account stress singularity near the 
crack tip. Variation of dynamic stress intensity factor with time and variation of 
dynamic fracture toughness with velocity were studied and compared with available 
experimental results. 

Nishioka and Atluri [2] gave elaborate information on the analysis of dynamic 
fracture using FEM. The most common ways to deal with the crack tip region was to 
simulate crack growth through gradual release of elemental nodal forces or 
embedding a moving element in which interpolation functions are determined by the 
continuum near tip fields at the crack tip in the mesh, and J-Integral consideration. 
Following special description the differential equation in time for the node point 
variables must be integrated. Because the dynamic fields associated with rapid crack 
growth are rich in high frequency content, small time steps are required for accuracy. 
Because of natural restriction to small time steps, the author’s report that many times 
combination of conditionally stable explicit time integration scheme and diagonalised 
mass matrix are found to be accurate. 

Chiang [6] presented a numerical procedure based on eigen function to 
determine the dynamic stress intensity factor for a crack moving at steady state under 
antiplane strain condition. An edge crack problem and a radial crack problem were 
solved using traction and displacement boundary conditions separately. It was shown 
that the dynamic effect was relatively insignificant for low crack propagation speeds 
provided that specimen size was fairly large. 

Thesken and Gudmundson [7] worked with elasto-dynamic moving element 
formulation incorporating a variable order singular element to enhance the local crack 
tip description. Kennedy and Kim [8] incorporated micropolar elasticity theory into a 
plane strain finite element formulation to analyse the dynamic response of the crack. 
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Materials with strong micropolar properties were found to have significantly lower 
dynamic energy release rate than their classical material counterparts. Wang and 
Williams [9] investigated high-speed crack growth in a thin double cantilever beam 
specimen using FEM. The cantilever end was loaded under an initial step 
displacement of 1mm and then pulled with a constant velocity the crack was 
propagated at assumed speed by gradual node release technique. Large dynamic 
effects were observed because of wave reflections within finite specimen size. The 
reflection and interaction of the waves from free boundary were avoided by limiting 
the duration of study. Beissel, Johnson and Popelar [10] presented an algorithm, 
which allows crack propagation in any direction but doesn’t require remeshing or the 
definition of new contact surfaces. This is achieved by tracking the path of the crack 
tip and failing the element crossed by the path such that they can no longer sustain 
tensile volumetric stress. The edges of these failed elements simulate crack faces that 
can sustain only compressive normal traction. Lin and Smith[ll] presented a method 
wherein crack growth is predicted in a step by step basis from Paris law using stress 
intensity factor calculated by using finite element method. The crack front is defined 
by a cubic spline curve from a set of nodes. Both the one quarter node crack opening 
and 3-D J-Integral method are used to calculate stress intensity factor. Automatic 
remeshing of finite element model to a new position that defines the new crack front 
enables the crack propagation to be followed. 


1.4 Present Work 

Present work is a modification of ‘Stiffness Release Model’ for dynamic crack 
propagation in DCB specimen( Siva Reddy [12] ). In the ‘Stiffness Release Model’, 
instead of reducing crack tip nodal force as is done in force release models, an 
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additional one dimensional elastic element is attached at the crack tip node whose 
stiffness value is gradually reduced, as the crack advance from one node to next node. 

In the present work, stiffness release model is further improved and the mass 
of the element surrounding the crack tip is modified as the crack propagates and the 
effect of various other parameters with crack velocity is also studied. 

Chapter 2 describes the basics of FEM formulation and the existing 
propagation model. 

Chapter 3 presents the crack propagation element model. 

Chapter 4 describes the results and discussion. 

Chapter 5 presents the conclusion of the present work and scope of future 

work. 
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CHAPTER 2 


BASICS OF FEM AND EXISTING CRACK 
PROPAGATION MODELS 

2.1 FEM Formulation 

2.1.1 FEM Equation 

Equations that govern the dynamic response of a structure are derived by using virtual 
work principle involving work of external forces, internal, inertial and viscous forces for 
any small admissible motion. 

The finite element solution involves discretization of domain (Q) into suitable 
elements and approximating the field variable interior to the element in terms of nodal 
values using suitable shape function. 

FEM equations can be derived from the virtual equation for the body as 


|{5s} r {a}dF+ \{5g} T p{q}dV+ \{5q) T k d \q'ftV (2. 1) 

q q n 

- {{5?} r {F,}rfF+ 

a r a 

where 

(&?!.{&( 

m) 

M 

p 

K 


: Virtual displacements and corresponding strains 
: Velocity and acceleration component 
: Stress component 
: Material density 
: Damping coefficient 
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\P B } : Body forces 
{P s } : Surface traction 
T CT : Surface traction boundary 


Displacement [q } (e) within element can be expressed in terms of its nodal values {(9} <e> as 

{q} (e) =[N] M {Qf e) (2.2) 

Strain-displacement relation can be written as 

{s} (e) =[5] (c) {0r (2.3) 

Stress-strain relation or constitutive equation is given as 

{cx} (e) = [Df * ' {s} (e) = [D] (e) [B] (e) {< 9} Ce) (2. 4) 

where, 

[iV] (e) - Matrix of shape function 
[s] (e) - Strain displacement matrix 
[ J D] <e) - Matrix of material properties 

The present analysis is made for an undamped system and neglecting body forces. 
The stresses, strains and displacements are expressed in terms of nodal variables as given 
in equations 2.2, 2.3 and 2.4. This on substitution in Eqn (2.1) gives, 


we prw'Wdv {«}+s -o 


(2.5) 


where summation is over all elements, N E . 
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since {5^} is taken arbitrary, satisfying kinematic conditions, the equation can be written 
as, 


If \p[N}' )T [N}y\d^dr\q} 

e=l V^Q(e) J e=\ y^e) J 



p}- >r {p.rvwn 


tF 


J 


( 2 . 6 ) 


The final set of assembled equations can be written as, 

M<2} + [M]jo}={F} (2.7) 

where 

[m]= ^[m] <£) , Global mass matrix 
[&"] = ^ [AT ] (e) , Global stiffiiess matrix 

{f} = , Global applied force vector 

and 

[iTj c) = [Dj e) \B J e) \j\dB,dr] 5 Elemental stiffiiess matrix 

n M 

[M| e) = | p[/V} e) [//J | j\dc,dr \ , Elemental mass matrix (2. 8) 

n“> 

^ = s' } \j\d^ , Elemental force vector 

j~(e) 


2.1.2 Integration algorithm 


Equations of equilibrium governing linear dynamic response of a system of finite 
element is given by Eqn (2.7). 


9 



Eqn 2.7 can be written as. 


F,(t) + F E {t) = Fit) (2.9) 

where, F, (t) are inertia forces, F,(t) = \M ]{q|, F e (t) are elastic forces , F E (t) = [X'jfg} , 
all of them being time dependent. Thus, in dynamic analysis, the displacements {<2},{e} 
and {F} are all functions of time. 

Mathematically, for small deformation A t can be assumed. Eqn (2.7) represent a system 
of linear differential equation of second order and the solution to the equation can be 
obtained by standard procedures of solving differential equations. 

The procedures for solving differential equations are divided into direct 
integration and mode superposition of which former is preferred in wave propagation 
problem. In this integration scheme, there are many different methods, which can be 
classified as ‘Explicit’ or ‘Implicit’ whose advantages and disadvantages are given by 
Bathe[13], 

In the present work two different integration schemes are used 

1 . N ewmark ’ s method 

2. 6 -Method 

which are discussed in the following Section. 

2.1.2.1 Newmark’s Method of Integration 

Newmark’s method is an extension of ‘Linear acceleration method* in which, linear 
variation of acceleration is assumed from time t to t + At . 


Velocity and acceleration approximated in terms of displacements are 

'“{9}='fel + l(l-«)'{«}+8""fe}W (2. 10) 
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( 2 . 11 ) 


t+At 


{qy {q}+‘ {q}At + 


(\ 


— a 


u 




AT 


where a and 8 are the parameters chosen suitably to have accuracy and stability. 
Usually, a and 5 are taken as 0.25 and 0.5 respectively to get unconditional stability. As 
mentioned earlier, equilibrium at t + At is considered along with these approximations. 
Solving equations 2.8 and 2.9, we get expressions for t+At {q] and t+Al {q}, in terms 
unknown displacement ,+iJ {q } . These are then substituted in the following equation to 
solve for <+A ' [q\ 

M'*" {$} + [£]'*" {?}='*" {F} (2. 12) 


Entire method can be summarised as follows. ( Bathe [13] ) 

1 . Form global stiffness matrix [K], and mass matrix [M], 

2. Initialise 0 {# }, 0 }, ° {?} . 

3. Select time step At and calculate integration constants 




A A 5 ^1 5 l - A "l 

a At 2 a At “ a At 2a 


; a 3 


< 2 4 — 1 , a 5 — 


At (5 


a 


2 \a 


— ~2 


; a 6 = At(l - 8 ) ; a n - 8At 


For effective stiffness, \K '] = [■&"]+ a Q [. M ] 


4. For each time step : 

(i) Calculate effective load at t + At 

,+Ar {f}+ [ Mja‘ 0 {q) + a[ {<?} + a[ fe}) 


(2. 13) 


(2. 14) 


(2. 15) 
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(ii) Solve for displacement at time t + At by [F'J +Ai {q]=‘ +At {F'} 

(iii) Calculate acceleration and velocity at time t + At 

‘ +A ‘ {4} = a Q (' +Af {q}- {9})- a‘ 2 {q\- a\ {?} 


t+&t 


{4}=' {q} + a[{q} + a 1 ;* {q) 


(2. 16) 


2. 1.2.2 0 -Method 


Hoff and Pahl[14] presented an unconditionally stable second-order-accurate 0 method. 
This method is found to be more stable compared to Newmark’s method as also will be 
shown later. 

0 - method can be summarised in brief as follows. 

1. Choose parameters 0,(O.95 to 1.0) and 0 O (1.0) 

2. Form global stiffness matrix [K], and mass matrix [M]. 

3. Initialise 0 {^ }, ° K ° {^' } • 

4. Select time step At and calculate integration constants 


4.0x0 2 (1.5 -0,)x 4.Ox0 2 

a o = — rr^ ; a \ = — K ; 

At At 

For effective stiffness, [AT'] -\k\+ a G \m\ 

5. For each time step : 

(iv) Calculate effective load at t + At 

M A/fe}+ 0.5 x At 2 {?})- {?}) 

+ (l.O-0,)A(fe}'-“ 

(v) Solve for displacement at time t + At by {q]= +ht {F'} 

(vi) Calculate acceleration and velocity at time t + At 


(2. 17) 

(2. 18) 


(2. 19) 
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,+A ' {?} = a o (' +A/ {<?}-' fe})~ X fe} + (l -0 - 20/ )fe} (2. 20) 

' +A< {gjs' fe}+ A t(Q l - 0.5)' fe} + (1.5- 0, )***" fe} 

0 -Method offers an unconditionally stable solution. While determining resultant force 
from the displacement boundary condition, by solving differential equation for one 
element, Newmark’s method was observed to be giving unstable oscillations in the 
solution. On the other hand, 0 -method was found to be giving stable solution under 
similar conditions. This is discussed in detail later in Chap 4. 

2.2 Methods for Simulation of Dynamic Crack Propagation 
by FEM 

To simulate a crack propagation in solids two different concepts of finite element 
modelling axe in use i.e., stationary mesh procedure and moving mesh procedure. Out of 
these two, moving mesh procedure, as the name implies involve change of mesh in each 
step as the crack propagates. This is computationally very intensive. Stationary mesh 
procedures do not require changes in the mesh in general, and require only change in the 
boundary and loading conditions. The two procedures are explained in the following 
Sections. 

2.2.1 Stationary Mesh Procedure: 

In the simple stationary mesh procedure of modelling linear elastodynamic crack 
propagation, the nodes ahead of the crack tip are spaced at V c At (V c being crack velocity) 
and crack propagation is simulated by releasing the node one at a time. However, if a 
simple node release technique is used, release of constraint on the displacement (at the 
preceding crack tip) one in one time step induces spurious high frequency oscillations in 
the finite element solution. To overcome these difficulties, several algorithms have been 
suggested in the literature to release the node gradually over few time steps. 
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Keegstra et al.[15,16] suggested a model in which node was not released 
instantaneously when the force at the crack tip reaches a critical value Fd, which is 
proportional to dynamic fracture toughness Kd. Instead the force was reduced from Fd 
gradually in accordance with an elastic spring model so that the nodal force vanish when 
critical displacement is reached. Another node release technique of gradual release of 
nodes is ‘Force release model’. 


2.2. 1.1 Force Release Model 

Suppose that the actual tip is located at ‘C’ in between the finite element nodes B and D 
as shown in Fig 2.1. The length segments BC and BD are b and d respectively. The 
holding back force ,F, at node B is gradually reduced to zero over a number of time steps 
as the crack tip reaches to the node D. Various schemes available to decrease the force to 
zero are as follows: 

(i) Malluck and King[17] suggested the release rate based on constant stress intensity 
factor 


F 


1 - 


( 2 . 21 ) 


where Fo is the original reaction force when the crack tip was located at node B. 

(ii) Rydhom et al.[l 8] Suggested the release rate based on constant energy release 


rate 


f 


,\ 


1 __ 

V dj 


( 2 . 22 ) 


(iii) Kobayasi et al.[19] suggested the linear release rate based on no physical 
argument other that pure intuition. 


F_ 

F n 






(2. 23) 


'■J 


In order to have more gradual and smooth propagation of crack Kishore, Kumar and 
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Verma [20] used a modified method. Holding back force at the crack tip B is linearly 
decreased to zero when the crack reaches end of the next element. Thus when crack tip 
goes beyond node B then 


F b 


H3 


f,_U 

d + d x ) 


(2. 24) 


where F HB is the reaction at node B, when the node was closed, b is the crack extension 
and d,di are the elements length as shown in Fig 2.1. And when the crack moves beyond 
node D to a point D) 


F e 


HB 


F n 


1 - 

V 

( 

1 - 


d + b 2 
d + d 


(2. 25) 


i J 


HD 


d x +d 


(2. 26) 


2 J 


where Fhd is the reaction at node D, when the node was closed- 

b,b 2 are the crack extension and d,di,d 2 are the element lengths as shown in Fig.2.1 


In force release model, time increment A t is usually taken such that it takes 15 to 
20 iterations for the crack tip to cross one element. The amount of force necessary to 
keep the nodes together is determined and a factor is used to proportionately decrease the 
force as the crack tip advances. Thus, at subsequent times, though the problem is highly 
dynamic, the current calculations were based on the force calculated to keep the nodes 
together 15 to 20 time steps before. This will cause discrepancy as the crack advances 
further from one element to the next one causing large oscillations in the solution. This 
discrepancy should be overcome for obtaining better solution. A new model is proposed 
to overcome the drawbacks of force release model which instead of considering dynamic 
force at the crack tip, takes into account the stiffness and mass of the element at the crack 
tip. It is described in detail in Chapter 3. 
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Fig 2.1 Crack opening scheme 
in force release model 


2.2.2 Moving Mesh Procedure: 


In moving mesh procedure, entire mesh moves with the crack tip. It can be further 
classified into an approach wherein a) the entire mesh moves b) only the mesh in a finite 
and small region around the crack-tip moves along with the crack-tip. Aberson et al. 
[21,22] used a singular crack-tip element as shown in Fig 2.2(a), which incorporates the 
first 13 terms of Williams’ eigenfunctions[23], appropriate to a stationary crack in a 
linear elastic body. The crack tip moves within the singular element, between the nodes A 
and B as shown in Fig 2.2(a). When the crack tip reaches the node B, the mesh pattern is 
changed suddenly, as illustrated in the figure. 

Another attempt at using Williams’ eigenfunctions was made by Patterson and 
01dale[24,25]. The singular element (Fig 2.2(b)) has 13 nodes and is topologically 
equivalent to two assembled 8 noded isoparametric elements. The location of the singular 
element is suddenly changed by a distance equal to the size of the regular element ahead 
of the crack tip, when the crack tip propagates a critical distance within the singular 
element. The displacements of the singular element are matched to those of the 
surrounding solid only at the nodes connected to the regular element. Thus, this model 
violates the displacement compatibility condition at the interfaces between the singular 
element and the surrounding regular elements. In another procedure, the singular element 
translates in each time step for which crack growth occurs. Thus crack tip is always at the 
centre of the singular element throughout the analysis. The regular elements surrounding 
the moving singular element are continuously distorted. Mesh pattern around the singular 
element is periodically readjusted for simulating large amount of crack propagation. 


17 



^ 

A' B' 


M- 




(a) (b) 

Fig 2.2 Discontinuously moving singular elements 
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2.3 Energy Release Rate and Its Determination 


In fracture mechanics a crack can be characterised by four different parameters 

i. Energy release rate ( G,,G n , G m ) is energy based and is.also applied to brittle and 
less ductile material 

ii. Stress intensity factor ( K , , K n , K m ) is stress based. This parameter is^alsO' applied 
to brittle and less ductile material. 

iii. J-Integral (J) which has been developed to deal with ductile material and can also be 
applied to brittle material as well. 

iv. Crack tip opening displacement (CTOD) is displacement based and is developed for 
ductile materials. 


From computational point of view, stress intensity factor require special element 5 
so as to simulate stress singularity at the crack tip whereas J-Integral and energy release 
do not require any such modelling of stress singularity. In energy release rate approach, 
energy is found out in entire domain and in J-Integral, evaluation is done along a path far 
away from crack tip. Thus both these approaches adroitly; avoid analysis close to the 
crack tip and don’t need any modelling of stress singularity. 


Griffith [26] outlined that in brittle fracture, the extension of a crack requires 
the formation of new surface, he r easoned with its associated surface energy, 
consequently, crack in a brittle solid should advance when the reduction of the total 
potential energy of the body during a small amount of crack advance equals the surface 
energy of the new surface thereby created. Most of the energy release, as the crack 


advances, comes from the parts of the plate which are adjacent to the cracked surface. 


- — Two important parameters need to be considered 

1) How much energy is released when a crack advance S 

2) Minimum energy required for crack to advance in forming two new surfaces 
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First parameter is measured in terms of energy release rate denoted by G. G is a 
function of crack size in general. Energy release rate is the amount of energy release per 
unit increase in area during crack growth. This energy is supplied by the elastic energy in 
the body and by the loading system. 

The energy requirement for a crack to grow per unit area extension is called as 
crack resistance and is usually denoted by symbol R. Crack resistance is sum of energy 
required to 1) form new surface and 2)cause anelastic deformation 
Both the available energy release rate and crack resistance are important to study the 
possibility of crack becoming critical. Obviously, when the available energy release rate 
far exceeds the crack resistance, the crack starts to grow at high speed. 

In the present study ‘energy release rate’ is adopted as it is a more comprehensive 
concept. 

As the crack advances, 

1) Stiffness of the component decreases. 

2) Strain energy in the component either increases or decreases. 

3) Work is done on the component if there is an application of load. 

4) Energy is being consumed to create two new surfaces. 

For quasi-static growth of crack by crack area AA, the incremental work by 
external force be , change in strain energy be AU, and the available energy, GAA, 
satisfy energy balance equation. 

GAA=AW ext -AU (2.27) 

G=~(£/-r„) (2.28) 


G=- 


dTI 

~dA 


(2. 29) 


where IT is the potential energy 

For a plate of uniform thickness, dA=Bda 



where da is the incremental crack length 
B is the thickness of the plate 


i an 

B da 


(2. 30) 


For dynamic crack propagation problem, the kinetic energy ,T, of the body should be 
taken into consideration. Dynamic energy release rate , Go, is different as some energy 
may be consumed to impart kinetic energy to the cracked portion of the body and to 
generate stress waves. 


For dynamic case, energy balance becomes 

GAA = AW ext -AU-AT (2.31) 

where AT is the increment in kinetic energy in the body. 


For constant velocity crack propagation, 

G = ~(W a -U-T) 

B da 

For crack moving with velocity v , 

Uw a -U-T) 

G = 

B v 

In the present study, Energy release rate is found out using equation (2. 17). 
Different energies are found out using 


(2. 32) 


(2. 33) 


= {Qf{F} 

(2. 34) 

\u]T-{qY[kIq} 

(2. 35) 

T T-{q} t [M ]{q} 

(2. 36) 
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where 

is the external work done 

U is the strain energy in the component 

T is the kinetic energy in the component 

{i 7 } is the external force vector 

[q\ is the displacement vector 

{9} is the velocity vector 

[Ai] is global stiffness matrix 

[M ] is global mass matrix 

2.4 Closure 

In this Chapter, formulation of FEM equation is described along with integration 
algorithms for solving equilibrium equation governing linear dynamic response of a 
system of finite element. Various methods of simulating dynamic crack propagation with 
stationary mesh and moving mesh are described. This Chapter also describes some 
general concepts in energy and energy release rate which will be used in the next 
Chapters to construct a propagation element. 
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CHAPTER 3 


CRACK PROPGATION ELEMENT 
3.1 Introduction 

Existing force release models have certain shortcomings, such as, 

1 . They give rise to very high oscillations in the solution. 

2. They are not suitable for solving inverse problem in crack propagation. If the 
crack propagates and stops in between an element and starts again later when 
another stress wave arrive, how much holding force to be used is not known. 

In the following sections, a propagation model with modified stiffness and mass 
matrices are presented to take into account of the crack extension. 


3.2 Stiffness Release Model 

In stiffness release model[12], one-dimensional elastic element is added at the crack 
tip node, whose stiffness is reduced gradually as crack tip advances to the next node 
(Fig 3.1). Additional stiffness required is very large when crack tip is at the beginning 
of element and zero when crack tip reaches to the other end. Amount of additional 
stiffness required is a function of crack length and is obtained from quasi-static crack 
propagation analysis. 

Fig 3.1(a) shows a symmetric part of deformed 2-D plate in which crack is at 
node B, gradually propagating to an intermediate position in between node B and 
node C 
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Let, 

K 0 be the stiffness at node B of original plate without additional spring 
u 0 be the displacement of node when crack tip reaches the next node and element 
opens up completely. 

K s can be found out from equilibrium equation as 


K 0 (u 0 - u ) = K s u 


(3. 1) 


K=K r 


— -1 

u 


Energy in the spring is given as 


E = —Ku 2 =-K t 

2 2 


(3. 2) 


(3.3) 


Energy release rate is written as 


G 


A E 
AA 


where A A is the increment in crack area 

Let a non-dimensional parameter for crack length (a) be defined as, 
a=crack length (within element) / d 
where d is the distance between two nodes. 

AA = BdAa 


(3.4) 


where B is the thickness of plate. 

For static case, energy release rate for a double cantilever beam can be written as 


G = 


1 

2 




1 du 
~Bd~da 


(3.5) 
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Assuming linear variation of G and a, G can be written as 

G = C, + C 2 a (3. 6) 

where, C, and C 2 are the constants to be determined. 


Substituting G in equation (3.5) 


C, 


+ C 2 & = — 


1 KqUq 

2 Bd 


(—) 

\daj 


(3.7) 


On integrating equation (3.7) 


2 

u = 

*0 «0 


Bd 


C x a + 




J 


On substituting 

initial condition u=0 when a=0 
end condition u=uo when a=l 


« o = “ 


* 0«0 


-Brf 


C, 


(3.8) 


(3*9) 


Eq.(3.8) and Eq.(3.9) gives 


h _ 2C,a + C 2 a 2 
Uq 2C, + C-y 


(3. 10) 



On simplification it reduces to 


where 


— = a + C(a -a 2 ) = / 
u n 



C 3 

2Cj + C 




2 J 


Eq. (3.2) finally can be written as 


K 


s 


= K o 


flzZl 

/ J 


( 3 . 11 ) 


( 3 . 12 ) 


( 3 . 13 ) 


where f=u/uo as given by Eq (3 . 1 1 ) 

Static stiffness K 0 can be determined for a given specimen under static condition. 
Energy release rate is determined using Eqn (3.5) and (3.11). 

However as Eqn 3.5 was based on quasi-static crack propagation, it doesn’t 
take into account the kinetic energy in the body. In dynamic problem total energy 
comprise strain energy in the body, kinetic energy and work done by external forces. 
For finding energy release rate, drop in the total energy in the system is considered in 
the present model and energy release rate is determined using fundamental Eqn (2.15) 

Effect of Constant ‘C’ 

Constant ‘C’ in equation (3.11) can be determined for quasi-static case so as to get a 
straight-line variation of energy release rate curve. For quasi-static case this value is 
very small. However it was observed that C has considerable effect in dynamic crack 
propagation. Effect of change in value of C on f Vs a curve is shown in fig (3.2) 
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Non Dimensional crack (a) 

Fig. 3.2 f Vs a curve for different C 


Spring stiffness ,K S , which is gradually reduced as the crack propagates, is a 
function of crack length ‘a‘ and constant C. This constant, C, is affected by change in 
element size and time step At. For a given problem, its value need to be chosen so that 
final energy release rate curve are smooth without wide fluctuation. 

Although the proposed stiffness K s of the additional elastic spring element 
varies from <x to 0, it changes discretely in a finite number of steps as the crack 
propagates through one element. For example, if the number steps to cross one 
element is considered as 10, then ‘a’ varies from 0.1 to 1.0 in steps of 0.1. Instead of 
finding K s at discrete points, an approach of averaging values over each step is 
adopted. This is brought about by averaging ‘f 4 in Eq. 3.1 1 in between for each small 
step and then finding out K s using Eq. 3.13. It is observed that averaging the values of 
f in each step gives smoother variation of energy release rate. This effect will be 
illustrated later in Chapter 4. 


28 





In order to bring the effect of mass variation, shape function Ni of the element at the 
crack tip is modified in such a way that it is a function of non dimensional crack 
length ‘a’ . 

In the present work modified Nj is proposed as, 




(1-0 

a 

( 1-0 

l 2 J 


v 2 j 


(3. 18) 


where 

Nj=l when^= 7 ] = -l 
=0 when E, = 77 = 1 


( ioq.i-O 0 - 9 
a J 


(3. 19) 


a is chosen such that, for a given crack length (a) within an element 
Ni=l when £= -1 
=0 when £>( 2 a-l) 
for example with a=0.5, 

Ni=l for £= -1 

=0 for £>1 or a>0.5 

Thus as the crack advances, shape function Nj varies in such a way that mass 
of the element increases gradually. Finally when a=l i.e., when crack is at the end of 
the element, shape function Ni take its bilinear form. Variation of N; along crack axis 
for various crack lengths is as shown in fig (3.5) and the surface plot of Ni for a=0.5 
is shown in fig (3.6). 
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Shape fun 
N, 


t 

Shape Fun 
N, 



Non dimensional crack length (a) 

Fig 3.5 Variation of shape function N] Vs 
crack length a (ri= -1) 



Fig 3.6 Shape fun Nj for a=0.5 




With modified shape function, element mass matrix given by Eqn (2.8) get 
modified. The effected coefficients are: 


q(c} V ^ / \ 2 / 


\j\d%dr\ 


(3. 20) 


M 12 = M n = 0 


(3.21) 


— M 3] — M 24 — M 42 — 


f! 

fi-0 

a 

ri-^ 


J 

(c) 

l 2 j 


l 2 J 

1 l 2 J 


J\dldr | (3. 22) 


M 14 = M 41 = M 23 = m 32 = 0 


(3. 23) 


m, 5 =m 5I =m 26 =m 62 = }( 

o(-)' 




V * J 


V 4 , 


\j\dqdri (3. 24) 


M 16 = M 61 = M is =M s2 =0 


(3. 25) 


A/"i 7 — M 71 — M 2g — M g2 — j" I 



a+l 


1,1 2 J 


V 4 J 


|/|^T] (3. 26) 


M 1S = M S1 = M 27 = M 11 = 0 


(3.27) 


Comparative figures of elemental mass matrices for a=0.5 and a=1.0 are shown Fig 
3.7 (a) and (b) respectively. 
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For a=1 .0 
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For a=0.5 
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(b) 


h^-8 Elemental mass matrix for a=1.0 and a=0.5 



By modifying shape function of elemental mass matrix, mass of the element is 
made to vary as a function of crack length. Element mass is zero when the crack is at 
the beginning of the node and increases gradually till crack reaches to the end of the 
element, when mass of the element becomes equal to the mass of a regular element. 
Effect of this modification on energy release rate variation will be shown in Chap 4. 


3.4 Closure 


Chapter 3 described in details of variation of various aspects of the proposed 
model. This model will be investigated with hypothetical data as well as experimental 
data in Chap 4. Effect of various parameters on energy release rate will be described 
in the succeeding chapter. 
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CHAPTER 4 


RESULTS AND DISCUSSION 
4.1 Introduction 

In this Chapter results of the dynamic crack propagation analysis are presented using the 
method described in the previous Chapters. First the method is validated for both quasi- 
static case and dynamic case. Later results are also presented for the available 
experimental data for crack propagation. 

Various aspects of the results involve, 

• Determination of C for quasi-static crack propagation. 

• Crack propagation at constant speed under dynamic input force pulse 

• Crack propagation at constant speed under time dependent input displacement (based 
on experimental data) 


4.2 Validation of the Proposed Model 


For the analysis purpose of investigating various aspects, a hypothetical problem 
of cantilever beam specimen is considered (Fig 4.1(a)). Due to symmetry only upper half 
portion is taken into account(Fig 4.1(b)). 

Material and geometric properties are as follows 

Modulus of Elasticity E=210xl0 9 N/m 2 

• Poisson ratio v=0.3 

• Length of specimen L=0.06m 
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• Width of specimen W=0.005m 

• Thickness of specimen B=0.024m 

• Crack length a'=0.03m 

The problem is symmetric and only one half is considered for descritisation and 
analysis. The DCB specimen is analysed by using 4-noded isoparametric element. Details 
are as follows 


• Number of elements 

600 

• Number of nodes 

671 

• Size of element in X-dir 

1.0mm 

• Size of element in Y-dir 

0.5mm 

For dynamic crack propagation. 

a short pulse is applied at the top of the specimen(Fig 

4.2). The pulse can be modelled by F(l-cos(cot)) whose details are as follows: 

• Pulse duration 

p=8ps 

• Force amplitude 

F 0 =50N 

• Frequency of pulse 

f=125KHz 


Time step At depends on the element size and wave velocity. 


4.3 Determination of C from Quasi-Static Case 


For quasi-static case, mass of matrix is essentially zero and it is only additional stiffness 
of elastic element, which will vary from maximum (infinity) to zero. Static stiffness K 0 

can be determined for a given mesh by applying force at the node preceding crack tip 
node and measuring corresponding displacement at that node. 
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For quasi-static case, a wedge loaded DCB specimen is considered with a 
constant displacement equal to 0.025xl0'frn at the crack opening. Energy release rate 
variation due to crack propagation is as shown in Fig 4.3 Curves are plotted for different 
values of constant C. It is observed that, for C= -0.01 and C=0.0 though G is continuous 
within the element, it is not continuous from element to element. On the other hand for 
C— -0.028, G curve is continuous both within the element and also across element 
boundaries. 

These C values may not be suitable for dynamic case. However, this forms the 
initial guess value. 


4.4 Dynamic Crack Propagation 

For dynamic crack propagation study, the same DCB specimen data as mentioned in Sec 
4.2 is considered. Crack propagation starts after the force pulse of 8ps is completely 
applied. 

For all the cases following parameters are taken 

i. Crack initiation time 8ps 

ii. Crack speed 1500m/s 

Two different time steps are considered. 

Before crack propagation Ati=l .Ox 1 0" 6 sec 

After crack propagation initiates At 2 a reduced value for higher accuracy 

4.4.1 Effect of Change in Value of C 

As obtained earlier in Sec 4.3 for quasi-static case, value of constant C for straight-line 
variation of energy release rate is C= -0.028. The same value when used for dynamic case 
gives energy release rate variation for crack velocity equal to 1500m/s as shown in 
Fig 4.4 


38 



Thus it can be seen from Fig 4.4 that there are large oscillations at the inter 
element boundaries. Besides this there are smaller oscillations within the element. It 
clearly shows that C value for the quasi-static case is not suitable for the dynamic case 
and has to be modified. Variation of C values causes effective stiffness value at the crack 
tip (in vertical direction), to vary as shown in Fig 4.5. C value is adjusted so that the 
energy release rate curve becomes smoother. Effect of change in C value on energy 
release rate curve is shown in Fig 4.6. It can be seen that C=0.4 gives rise to less 
oscillations in Gi in comparison to other values and it is chosen for further analysis. 

4.4.2 Dynamic Crack Propagation with Averaging of ‘f 

As illustrated in Sec 3. S. since only finite number of steps are considered for the crack 
propagation in one element. The zone between zero and first step goes unconsidered. So a 
method of averaging the values of ‘f in Eqn. 3.11 is considered. 

Fig 4.7 shows difference in energy release rate curve for the same value of constant C. 
Figure clearly shows that energy release rate curve with averaged value of ‘f ‘ gives 
better solution. 


4.4.3 Effect of Shape Function modification on Element Mass Matrix 

As the stiffness of elastic element decreases gradually to zero, the element mass can be 
increased from 0 to its usual value by modifying shape function Ni for element mass 
matrix for an element containing propagating crack (Sec 3.5). Energy release rate after 
this modification is as shown in Fig 4.8. Fig reveals that final energy release rate curve 
gives smooth variation within element and fairly stable solution across element 
boundaries. 

Similar results for crack velocity of 1250m/s are shown in Fig 4.9. It can be seen very 
clearly that solution with mass modification is better than one without mass modification. 
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4.4.4 Effect of element size and time step 


In the present model, for simulating dynamic crack propagation, both stiffness as well as 
mass matrix are modified for the element containing crack tip. For the crack tip node the 
effective combined stiffness and mass matrix value in the time integration scheme 
(Newmark’s method) is 

K' = K + ^—y-K 0 +a 0 M (4.1) 

= K + K s +a 0 M (4.2) 

where, 

K s is stiffness of elastic element and 
a 0 is a Newmark’s integration constant 


1 


a o = 


a At 


(4. 3) 


Thus effective stiffness K' depends on, K which in turn depends on size of 
element, K s which is a function of crack length and static stiffness value K 0 . And a 0 


depends on time step At , M which is a function of crack length ‘a’. Fig (4.10) shows the 
nature of (K + K s ) Vs non-dimensional crack length. ( K + K s ) is a monotonically 
decreasing curve and is independent of change in time step, At . Fig (4.1 1) show variation 
of a 0 M with ‘a’. This curve has increasing variation with crack length. One important 
point to be recognised is that decrease in At by half increases a 0 M four folds. Fig 4.12 

shows the combination of all terms giving by Eqn(4.2) and represents effective stiffness 
matrix at the crack tip node. 


Following observation can be made 

1. For At = 0.13/xs effective stiffness matrix has a monotonically decreasing nature 
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2. For At = 0.067 jjs , the curve decreases and then increases in the later half of the but 
otherwise it’s a smoothly decreasing curve. 

3. For At = 0.033fjs inertial term becomes dominant and effective stiffness is no longer 
a continuously decreasing curve. There is large sag in the middle of the curve. 

Fig 4.13 reveals its effect on energy release rate curve. 

1. For At = 0.13/15 because of higher time step inertial term has a negligible 
contribution and energy release rate has variation close to the case without mass 
modification case. 

2. For At = 0.067 /is there is a balanced contribution from mass and stiffness. Energy 
release rate curve gives smooth variation throughout crack length. 

3. For At = 0.033 /is oscillations are appearing both within the element as well as across 
inter element boundaries. 

These leads to'conclusion that size of the element and time step (At) should be chosen in 

such that K' decrease fairly monotonically. 


4.5 Crack Propagation in Glass Fabric/Epoxy DCB 
Specimen 

The present method is applied to determine the dynamic Gi for the experimental data as 
reported by Chowdhary[27] . 

From the mechanics point of view fibre composites are among the class of 
material called ortho tropic materials. Generalised Hook’s Law for 2-D ortho tropic 
laminate is given by Agarwal[28]. 

°ij=Qij s ij ( 4 - 4 > 

ij=l,2....,6 
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where cr ; axe the stress components, Q i} is the stiffness matrix, s y . are engineering strain 
components. 

For 2-D orthotropic material, 


where, 


q 


re.» 

Qll 

[ 

0 


f \ 

1 2 

► = 

021 

Q 22 

0 

< 

£ 2 > 

1512, 


_ 0 

0 

Q&6 _ 


y 12 . 


(4.5) 


fi,.= 


El 


1 v LT v 7 - L 


(4.6) 


e 22 = 


1 v LT v TL 


(4.7) 


Q\2 ~ 


^TL^L 

1 — ^lt^tl 


^_LT^T_ 

1 — ^ LT^TL 


(4.8) 


0-66 ~ G L T 

E l and E r are elastic modulae in Longitudinal and Transverse direction 
G lt is the shear modulus associated with axes of symmetry 
v LT and v TL are major and minor Poisson ratios. 


Relation between four of the five elastic constants is 


El 

E t 


^LT 

^TL 


Following data is used as input for finding energy release rate. 

• Displacement at cantilever end as a function of time 

• Velocity of crack propagation 

• Material data and specimen geometry 


(4. 9) 


(4. 10) 


42 



Material and geometric properties are as follows 


Material type 
Material density 
Elasticity' - constants 


Poisson ratio 


Glass fabric/epoxy laminate 

p=T825Kg/m 3 

El=26x10 9 N/m 2 

Er=6.6xl0 9 N/m 2 

GLT^S.SxlO^/m 2 

vlt^O.21 

v tl =0.05331 


3 different experimental data is used. 

Details of geometric data is given in Table 4. 1 

Table 4.1 Details of specimen 


Expt No 

L(mm) 

W(mm) 

B(mm) 

a' (mm) 

1 

7.8 

4.0 

24.7 

40.4 

2 

7.6 

3.8 

25.1 

40.3 

3 

7.6 

3.8 

24.8 

40.0 


L - Length of the DCB specimen 
W - Width of specimen 
B - Thickness of specimen 
a' - Crack length 

In the experiment, time history of displacement at the cantilever end is known. 
Fig 4.14, Fig 4.15 and Fig 4.16 gives displacement applied at the cantilever end as a 
function of time. For determining work done, it is essential to obtain equivalent force as a 
function of time applied at the cantilever end. Nodal force is obtained from finite element 
equation. 
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Resultant of equilibrium force obtained for a typical case is shown in Fig 4.17. It 
is evident that there are very high oscillations in the solution due to instability of the 
integration algorithm. 

A more stable scheme for solving differential equation the Theta method 
summarised in Sec 2.13 is used. The difference in the results by these two schemes are 
shown in Figs 4.17 and 4.18. Thus, oscillations in the equivalent force have reduced 
drastically. Small oscillations that are still left are finally eliminated by using smooth 
Bezier curve, which gives equivalent load curve as shown in Fig 4.19. 

Fig 4.20-4.22 give the dynamic Gi for these experimental data and results are 
summarised in Table 4.2 

Table 4.2 Initiation and propagation toughness with ‘Force release model’ and 
‘present’ model: 


Expt No 

V (m/s) 

T ini (s) 

G,„i (J/m 2 ) 

Gprop (J/m 2 ) 

FRM[27] 

Present 

FRM[27] 

Present 

1 

1233 

54.04 

418 

270 

11 

4.5 

2 

1367 

45.97 

116 

52 

9 

3.4 

3 

860 

44.86 

52 

19 

18 

5 


Energy release rate variation for experiment no 1,2 and 3 are shown in Fig 4.20, Fig 4.21 
and Fig 4.22 respectively. Similar results for force release model is shown in Fig 4.23, 
Fig 4.24 and Fig 4.25 

It is worth noting that energy release rate curve obtained by present analysis is 
much smoother and can be used for finding propagation toughness. 


44 






























Ffl-cosfcotY) 



Fn-cosfcotY) 


Fig 4.1 DCB Specimen (Full) 
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Fig 4.3 Energy release rate variation for quasi- 
static crack propagation 
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Fig 4.5 Effect of change in C on 
effective stiffness value at the crack tip 
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Fig 4.6 Effect of change in C on 
energy release rate 
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Fig 4.7 Effect of averaging ‘f values on 
energy release rate 
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Fig 4.9 Energy release rate for crack 
velocity 1250m/s 
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Fig 4.10 Effective stiffness at crack tip Vs ‘a’ 

without shane function modification 
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Non dimensional crack length (a) 


Fig 4.11 Variation of effective mass at crack 
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Fig 4.12 Variation of effective stiffness at 
crack tip Vs ‘a’ 
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Fig 4.13 Energy release rate for different 
time steps 
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CHAPTER 5 


CONCLUSIONS AND SCOPE FOR FUTURE WORK 

5.1 Conclusions 

A new model is proposed to investigate high speed crack propagation. This model 
involve further modification of stiffness release model and based on the results and 
discussion in Chap 4, following conclusions can be drawn. 

1 . Model gives fairly smooth and more stable variation of energy release rate. 

2. Experimental results reveal that unlike force release model energy release rate doesn’t 

drop to a low value in very few time steps but falls gradually which is a more 
practical result. 

3. Out of two methods for solving equation of equilibrium (Newmark’s and 
6 - method), 6 - method was observed to be more stable. 

5.2 Scope for Future Work 

1 . The crack propagation model may be further modified to get still better results. 

2. Present method can be extended for Mode II and mixed mode crack propagation. 

3. This method can also be extended for 3-Dimensional crack propagation problems. 
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